#include "fermiqcd.h"

void generate_stochastic(fermi_field& phi,
			 gauge_field& U,
			 coefficients& quark,
			 float absolute_precision,
			 float relative_precision) {
  int nc=phi.nc;
  fermi_field dummy(phi.lattice(),nc);
  set_random(dummy);
  mul_invQ(phi,dummy,U,quark,absolute_precision,relative_precision);
}


int main(int argc, char** argv) {
  mdp.open_wormholes(argc,argv);
  define_base_matrices("FERMILAB");

  char filename_gauge[128]="cold";
  char filename_n[128];
  char filename_m[128];
  int nt=16,ns=4;
  int nstochastic=10;
  coefficients light_quark;
  light_quark["kappa_t"]=0.11;
  light_quark["kappa_s"]=0.11;
  // light_quark["c_{sw}"]=0.1;
  float ap=1e-6;
  float rp=1e-6;
  default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;
#ifdef SSE2
  default_fermi_action=FermiCloverActionSSE2::mul_Q;
#else
  default_fermi_action=FermiCloverActionFast::mul_Q;
#endif
  int box[]={nt,ns,ns,ns};
  mdp_lattice lattice(4,box);
  gauge_field U(lattice,3);
  fermi_field source_n(lattice,3);
  fermi_field source_m(lattice,3);
  fermi_field sink_n(lattice,3);
  fermi_field sink_m(lattice,3);
  mdp_site x(lattice);
  mdp_site y(lattice);
  int i,j,k,n,m,a,b,t;
  mdp_complex sink;
  float pion[100];

  set_cold(U);
  for(n=0; n<nstochastic; n++) {
    generate_stochastic(source_n,U,light_quark,ap,rp);
    sprintf(filename_n,"%s_stochastic_source_%.3i", filename_gauge,n);
    source_n.save(filename_n);
    mul_Q(sink_n,source_n,U,light_quark);
    forallsites(x)
      for(a=0; a<4; a++)
	for(i=0; i<3; i++)
	  sink_n(x,a,i)=conj(sink_n(x,a,i));
    sprintf(filename_n,"%s_stochastic_sink_%.3i", filename_gauge,n);
    sink_n.save(filename_n);
  }

  float tmp=0;
  for(t=0;t<nt;t++) pion[t]=0;
  for(n=0; n<nstochastic; n++) {      
    cout << "n=" << n << endl;
      // read stochastic prop n
    sprintf(filename_n,"%s_stochastic_source_%.3i", filename_gauge,n);
    source_n.load(filename_n);
    sprintf(filename_n,"%s_stochastic_sink_%.3i", filename_gauge,n);
    sink_n.load(filename_n);

    forallsites(x)				
      for(a=0; a<4; a++)
	for(i=0; i<3; i++)
	  tmp+=real(source_n(x,a,i)*conj(source_n(x,a,i)));
    cout << "norm=" << tmp/(n+1) << endl;

    for(m=0; m<nstochastic; m++) {
      sprintf(filename_m,"%s_stochastic_source_%.3i", filename_gauge,m);
      source_m.load(filename_m);
      sprintf(filename_m,"%s_stochastic_sink_%.3i", filename_gauge,m);
      sink_m.load(filename_m);
      
      y.set(0);
      sink=0;
      for(a=0;a<4;a++)
	for(i=0;i<3;i++)
	  sink=sink_n(y,a,i)*conj(sink_m(y,a,i));

      forallsites(x) {
	t=x(0);
	for(a=0;a<4;a++)
	  for(i=0;i<3;i++)
	    pion[t]+=real((source_n(x,a,i)*conj(source_m(x,a,i)))*sink);
      }
    }    
  }
  for(t=0;t<nt;t++)
    cout << t << ", " << pion[t] << endl;

  mdp.close_wormholes();
}
  
